* --------------------- *
* This file loads the child level data to compute 
* the multiple-hypothesis p-values reported in Tables C1, C2 and C3.
* Author: Victor Ronda
* Date Created: 01 Mar 2020
* Last Update: 22 Jan 2023
* --------------------- *


clear all
set more off

**************************
*** Health and Dominance Measures ***
**************************

use "offspring_data.dta", replace

rename health_issue_*_10m3y hi_*_13
rename dominance* dom*
rename dom_relative_rank dom_pct


replace hi_m_13=1-hi_m_13
replace hi_m_13=. if obs_10m3y<=5

xtile mcoh1=mom_cohort if hi_m_13 ~=.,n(20)
xtile coh1=cohort if hi_m_13 ~=.,n(20)

local controls1 "coh1 female primip obs_10m3y"
local controls2 "mcoh1 mom_primip coh1 female primip"
local controls3 "mcoh1 mom_primip coh1 female primip obs_10m3y"

xtile mcoh3=mom_cohort if dom_pct ~=.,n(20)
xtile coh3=cohort if dom_pct ~=.,n(20)

local controls5 "coh3 female primip agedom"
local controls6 "mcoh3 mom_primip coh3 female primip"
local controls7 "mcoh3 mom_primip coh3 female primip agedom"

gen subgroup=1 if rear_inte==0|rear_inte==1
replace subgroup=2 if rear_inte==0|rear_inte==2
replace subgroup=3 if rear_inte==0|rear_inte==3
replace subgroup=4 if rear_inte==1|rear_inte==2
replace subgroup=5 if rear_inte==1|rear_inte==3
replace subgroup=6 if rear_inte==2|rear_inte==3

keep if rear_inte~=.

gen treat0=(rear_inte==0)
gen treat1=(rear_inte==1)
gen treat2=(rear_inte==2)
gen treat3=(rear_inte==3)

local controls ""
local controls2 ""

********************
***** Table C1 *****
********************

*** MHT(1) ***
* - The values in the third column (Romano-Wolf p-value) are the p-values for MHT(1)
* - The values in the second column (Resample p-value) are the Bootstrap p-values.

rwolf2 (reg hi_m_13 treat3 treat0 treat2 `controls')  (reg dom_pct treat3 treat0 treat2 `controls2') , indepvars(treat3, treat3) seed(3) reps(500) cluster(mom_id)
rwolf2 (reg hi_m_13 treat2 treat1 treat3 `controls')  (reg dom_pct treat2 treat1 treat3 `controls2')  , indepvars(treat2, treat2) seed(3) reps(500) cluster(mom_id)
rwolf2 (reg hi_m_13 treat3 treat0 treat1 `controls')  (reg dom_pct treat3 treat0 treat1 `controls2')   , indepvars(treat3, treat3) seed(3) reps(500) cluster(mom_id)
rwolf2 (reg hi_m_13 treat1 treat2 treat3 `controls')  (reg dom_pct treat1 treat2 treat3 `controls2')  , indepvars(treat1, treat1) seed(3) reps(500) cluster(mom_id)
rwolf2 (reg hi_m_13 treat2 treat0 treat3 `controls')  (reg dom_pct treat2 treat0 treat3 `controls2')  , indepvars(treat2, treat2) seed(3) reps(500) cluster(mom_id)
rwolf2 (reg hi_m_13 treat3 treat1 treat2 `controls')  (reg dom_pct treat3 treat1 treat2 `controls2') , indepvars(treat3, treat3) seed(3) reps(500) cluster(mom_id)



*** MHT(2) ***
* - The values in the third column (Romano-Wolf p-value) are the p-values for MHT(2)
rwolf2  (reg hi_m_13  treat3 treat0 treat2  `controls') (reg hi_m_13 treat2 treat1 treat3 `controls') (reg hi_m_13 treat3 treat0 treat1 `controls')  (reg hi_m_13 treat1 treat2 treat3 `controls') (reg hi_m_13 treat2 treat0 treat3 `controls') (reg hi_m_13 treat3 treat1 treat2 `controls'), indepvars( treat3, treat2, treat3, treat1, treat2, treat3) seed(3) reps(500) cluster(mom_id)
rwolf2(reg dom_pct  treat3 treat0 treat2  `controls2') (reg dom_pct treat2 treat1 treat3 `controls2') (reg dom_pct treat3 treat0 treat1 `controls2')  (reg dom_pct treat1 treat2 treat3 `controls2') (reg dom_pct treat2 treat0 treat3 `controls2') (reg dom_pct treat3 treat1 treat2 `controls2'), indepvars( treat3, treat2, treat3, treat1, treat2, treat3) seed(3) reps(500) cluster(mom_id)


*** MHT(3) ***
* - The values in the third column (Romano-Wolf p-value) are the p-values for MHT(3)
rwolf2 (reg hi_m_13  treat3 treat0 treat2  `controls') (reg hi_m_13 treat2 treat1 treat3 `controls') (reg hi_m_13 treat3 treat0 treat1 `controls')  (reg hi_m_13 treat1 treat2 treat3 `controls') (reg hi_m_13 treat2 treat0 treat3 `controls') (reg hi_m_13 treat3 treat1 treat2 `controls')  (reg dom_pct  treat3 treat0 treat2  `controls2') (reg dom_pct treat2 treat1 treat3 `controls2') (reg dom_pct treat3 treat0 treat1 `controls2')  (reg dom_pct treat1 treat2 treat3 `controls2') (reg dom_pct treat2 treat0 treat3 `controls2') (reg dom_pct treat3 treat1 treat2 `controls2'), indepvars(treat3, treat2, treat3, treat1, treat2, treat3, treat3, treat2, treat3, treat1, treat2, treat3) seed(3) reps(500) cluster(mom_id)
 
 
 
**************************
*** Fertility Measures ***
**************************
use "mother_data.dta", replace
keep if female==1

gen primip=0 if pregorder~=.
replace primip=1 if pregorder==1
egen numpreg_p1=rowmax(pregorder_ch1 pregorder_ch2 pregorder_ch3 pregorder_ch4 pregorder_ch5 pregorder_ch6 pregorder_ch7 pregorder_ch8 pregorder_ch9 pregorder_ch10 pregorder_ch11 pregorder_ch12 pregorder_ch13 pregorder_ch14 pregorder_ch15 pregorder_ch16)
gen numpreg=numpreg_p1
replace numpreg=0 if numpreg==. & female==1
replace numpreg=. if female==0
gen anyoffspring=1 if numpreg>0 & numpreg~=.
replace anyoffspring=0 if numpreg==0
gen age1preg2=(dob_ch1-dob)/365


gen anyoffspring_a8=anyoffspring if age2017>=8
gen numpreg_a8=numpreg if age2017>=8
gen numpreg_p1_a8=numpreg_p1 if age2017>=8 

local controls "cohort primip"

gen MR_vs_SPR=1 if rear=="MR"
replace MR_vs_SPR=0 if rear=="PR"|rear=="SPR"

drop if cohort==. | female==. | primip==. | MR_vs_SPR==.


label var MR_vs_SPR "$\Psi(s_{*,1})-\Psi(s_{*,0})$" 

********************
***** Table C2 *****
********************

*** MHT(1) ***
* - The values in the third column (Romano-Wolf p-value) are the p-values for MHT(1)
* - The values in the second column (Resample p-value) are the Bootstrap p-values.

rwolf2 (reg anyoffspring_a8  MR_vs_SPR `controls') (reg age1preg  MR_vs_SPR `controls') (reg numpreg_a8  MR_vs_SPR `controls') (reg numpreg_p1_a8  MR_vs_SPR `controls') , indepvars(MR_vs_SPR, MR_vs_SPR, MR_vs_SPR, MR_vs_SPR) seed(3) reps(500) 



*******************************
*** Birth Outcomes ***
*******************************
use "offspring_data.dta", replace

replace stillborn=0 if stillborn==.

foreach var of varlist bw diedb1m {
replace `var'=. if stillborn==1
}

gen livebirth=1-stillborn
gen liveda1m=1-diedb1m

local controls1 "cohort female primip"
local controls2 "mom_cohort mom_primip"
local controls3 "cohort female primip"

gen mom_MR_vs_SPR=1 if mom_rear=="MR"
replace mom_MR_vs_SPR=0 if mom_rear=="PR"|mom_rear=="SPR"

label var mom_MR_vs_SPR "$\Psi(s_{1,*})-\Psi(s_{0,*})$" 

drop if mom_cohort==. | female==. | mom_primip==. | cohort==. | primip==. |mom_MR_vs_SPR==.

********************
***** Table C3 *****
********************

*** MHT(1) ***
* - The values in the third column (Romano-Wolf p-value) are the p-values for MHT(1)
* - The values in the second column (Resample p-value) are the Bootstrap p-values.

rwolf2 (reg livebirth mom_MR_vs_SPR mom_cohort mom_primip cohort female primip) (reg liveda1m mom_MR_vs_SPR mom_cohort mom_primip cohort female primip) (reg bw mom_MR_vs_SPR mom_cohort mom_primip cohort female primip) , indepvars(mom_MR_vs_SPR, mom_MR_vs_SPR, mom_MR_vs_SPR) seed(3) reps(500) cluster(mom_id) 



 